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Abstract. We derive a purely algebraic framework for the identification of hierarchy 
equations of motion that induce completely positive dynamics and demonstrate the 
applicability of our approach with several examples. We find bounds on the violation 
of complete positivity for microscopically derived hierarchy equations of motion and 
construct well-behaved phenomenological models with strongly non-Markovian revivals 
of quantum coherence. 
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1. Introduction 

Every quantum system inevitably interacts with its environment. This interaction 
contributes to the emergence of classicality, represents the major challenge in the 
development of quantum information technological hardware [1], but can also be 
exploited for advanced means to control quantum systems like laser cooling [2] or 
dissipative state preparation [3]. 

The central difficulty in describing open quantum systems resides in the mere 
size of an environment. Most approaches aim at the reduction of such environments 
to their aspects that are most relevant for the system dynamics. This ranges from 
effective descriptions in terms of the system’s degrees of freedom only [4, 5], via 
approaches that explicitly include the most relevant environmental degrees of freedom 
[6, 7, 8] to numerically expensive treatments that aim at a potentially exact description 
[9, 10, 11, 12, 13]. In practice, one needs to take a compromise between accuracy and 
numerical effort; a comparatively simple effective description can help to develop a more 
intuitive model, but may also fail to identify more subtle features. 

A completely different approach is based on a phenomenological modelling of open 
quantum dynamics that promises a better understanding and deeper insights into the 
underlying physics. The hallmark is provided by the Kossakowski-Lindblad equation 
[4, 5] that allows us to construct models that respect the probabilistic interpretation 
of quantum mechanics: Markovian Lindblad master equations ensure the property of 
complete positivity such that all quantities that are interpreted as probabilities are 
always non-negative. Furthermore, it permits an interpretation of the environmental 
effects in terms of rates and elementary processes. The intuitive comprehension that 
comes along with this model suffers the disadvantage of a strongly limited applicability. 
Many systems that are actively investigated, like quantum dots [14], light-harvesting 
systems [15], or photonic band-gap materials [16], are characterized by a rather strong 
system-environment coupling giving rise to non-Markovian effects like the back-flow of 
information that are not captured by this model. 

For that reason, it is desirable to find a framework for the description of non- 
Markovian processes that ensures complete positivity and is in the best case even 
phenomenologically accessible. Complete positivity is not only crucial for physically 
sensible dynamics, but understanding the very conditions under which the equation of 
motion induces completely positive (CP) processes must also be considered a crucial step 
towards a phenomenological understanding. However, the addition of memory to the 
system dynamics gives rise to substantial complications [17] in particular with respect 
to complete positivity. Most approaches to this problem focus on integral equations over 
a memory kernel and for certain classes like semi-Markov processes [18, 19], collision- 
model-based approaches [20] or continuous time quantum random walks [21] conditions 
for complete positivity have been found. Some of these classes can also be covered by 
generalized sufficient conditions that have been obtained recently [22]. Eventually, the 
post-Markovian master equation [23] is capable to describe CP non-Markovian dynamics 
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by interpolating between the generalized measurement interpretation of the Kraus 
operator sum and the notion of a continuous measurement for Markovian processes 
but eventually also imposes conditions on the memory kernel. 

In this article we propose a different take on the question of CP non-Markovian 
dynamics. Rather than by integral forms over a memory kernel, our approach is 
inspired by the framework of hierarchical equations of motion (HEOM) [9] that employ 
- in addition to the system density matrix gi{t) - a set of auxiliary operators gi{t) 
{i = 2,... ,n) that contain information about the environment and system-environment 
correlation. Equations of motion of the form 


n 



( 1 ) 


j=0 


have been derived based on microscopic models for a wide range of systems including 
light-harvesting complexes [24, 25], molecular transistors [26], as well as electrons in 
tunneling junctions [27] and organic heterojunctions exposed to a phonon bath [28]. 
In particular they have proven to be very reliable in the regime where comparable 
intra-system and system-environment interactions prevent approximations based on 
separations of time-scales. 

We will start out assuming the structure (1) with a finite-dimensional system state 
gi(t) and target the identification of conditions on Cij such that the dynamical map Ai(t) 
with t > 0 defined via gi{t) = Ai(f) pi(0) is CP. In the case of n = 1, i.e. equations of 
motion for the system state only, it is well established that this condition is satisfied if 
and only if Cn can be written in Lindblad form [4, 5]. A generalization of this criterion 
to hierarchy equations with n > 1 that also permits an incorporation of non-Markovian 
dynamics is, however, not known. 

In this article, we propose a new method by means of which sufficient conditions 
for CP dynamics can be derived for HEOMs. In some cases, such conditions can be 
obtained analytically, whereas highly efficient numerical methods are in place whenever 
an analytical solution is not feasible. In section 2.1, the formal and algebraic framework 
is set before a very intuitive geometrical interpretation of the procedure is given in 
section 2.2. In section 2.3, we provide practical tools that ease an analytical examination 
as it is carried out for a couple of examples in section 2.4. More advanced techniques 
that aim for a numerical investigation of more complex problems are subject to section 3, 
where we introduce a transformation that permits a universal applicability of the method 
(section 3.1) and demonstrate how semi-definite programming can be employed for a 
highly efficient numerical examination of HEOMs (section 3.2). A truncation of a HEOM 
is often carried out in order to obtain approximate dynamics but can also give rise to 
equations that violate complete positivity. That is why in section 4.1 and 4.2, our 
approach is modified such that the binary classification into CP and non-CP dynamics 
is replaced by a quantification of the violation of complete positivity. This permits the 
application of our method to an important example of a truncated HEOM that describes 
the spin-Boson model in section 4.3, before we conclude our results in section 5. 
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Since complete positivity is a property of the dynamical map rather than of the system 
state, it is helpful to re-formulate (1) in terms of time-dependent maps Aj(t) that satisfy 
Qi{t) = Ai{t) Pi(0). The equations of motion for these dynamical maps are obtained from 
(1) and read 

n 

^jit) , (^ = 1, 2,..., n) . (2) 

This set of equations can be expressed in a very concise form, when we introduce the 
extended dynamical map A = ® N) with {|i)} being the orthogonal and 

normalised basis of an n-dimensional Hilbert space. This notion of extended objects 
that are formed by operators is highly convenient and will in the course of this article 
be denoted by bold-faced symbols. The scalar product on this extended space reads 

{A,r)5tr(A'r) = ^tr(Alr,) (3) 

i 

for any two extended vectors A and F. In the same spirit, we also introduce the extended 
state git) = ® |i) and dehne the extended generator C = Eij 0 \i){j\, 

which permits to express (2) in the compact form 

dtA{t)=CA{t) . (4) 

By means of this notation we are in the position to concisely describe the algebraic 
framework for the derivation of conditions for completely positive dynamics. In practice, 
it is often convenient to work with Bloch-type representations, where the Bloch vector 
® N) of 8 is dehned in terms of the Bloch vectors bi of the operators Qi, i.e. vectors 
with components [bi]j = tr^ajQi) for j = 0,x,y,z, where ax, ay and cr^ are the Pauli 
matrices and cxo = I, or their generalisation to higher dimensional systems. For extended 
maps A and generators jC the corresponding representations are dehned analogously. 

2.1. The algebraic framework 

Let us hrst detail the concept of valid maps. An extended dynamical map A(t) is valid 
if and only if the system map Ai(t) is CP, i.e. it can be written as 

Ai(t) = X] Xij{t) AiQiA] , (5) 

ij 

with a positive semi-dehnite matrix x{t) and a set of mutually orthonormal operators 
fii. Showing complete positivity of Ai(f) is thus equivalent to proving positive semi- 
dehniteness of x(t). 

The initial map Ai(0) at the time t = 0 is given by the identity map such that 
x(0) has rank one. The rank of a matrix provides information about the number of 
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non-vanishing eigenvalues and does therefore also indicate the order of the highest non¬ 
vanishing elementary symmetric polynomial. Those elementary symmetric polynomials 
read 

ei(x) = tr(x), e 2 (x) = [tr^(x) - tr(x^)]/2, • • ■, e^(x) = det(x) , (6) 

where r x r is the dimension of y, and when y has rank k, then ej{x) = 0 for all J > k. 
For the initial matrix x(0) one obtains ei(x(0)) = 1 and ej(x(0)) = 0 for all j > 1 but as 
time evolves, the eigenvalues of x{t) will typically take non-vanishing values such that 
the rank of x{t) increases. We will now consider a time fp > 0 at which all non-trivial {i.e. 
not permanently vanishing) eigenvalues of x{t) Eire strictly positive. This is equivalent 
to x{tp) being positive semi-definite and having rank h = max(rank(x(f)) 11 > 0). 
Any sign-flip of any eigenvalue of x(f) for t > tp is then necessarily accompanied by 
a vanishing polynomial Ch, for why it is sufficient to ensure that Chixif)) does never 
become equal to zero for t > tp in order to show complete positivity from time tp 
on. If no eigenvalue of x(f) remains equal to zero {i.e. h = r), then Ch is given by 
the determinant but even when one or more eigenvalues of x(f) vanish permanently, 
the proposed method will still be applicable with Ch being an elementary symmetric 
polynomial with lower degree. 

In order to prove eh{x{t)) > 0 for t > tp, it is assumed (the assumption will be 
justihed in section 3.1) that one can hnd a non-linear transformation of the coordinates 
A{t) such that eh{x{t)) is quadratic in the extended map A(t) and can be written as 

en{x{t)) = (A(0),5A(0)) - (A(t),5A(t)) (7) 

with a suitable operator S and satisfy S' = Sb 

As we will show in the following the dynamical map Ai(f) is CP for all times t > tp, 
if there exists an operator R = Rl that 

(i) fulhls the operator inequality R + RC < 0, 

(ii) is “normalized” as {\{tp), {R — S) A{tp)) = 0, 

(iii) is such that i? — S is positive semi-definite. 

To prove this statement, we employ (ii) to reformulate (7) to 

ehixit)) = ehixiip)) + (A(tp), A(tp)) - (A(t), S A(t)) . (8) 

Since eh{x{tp)) is strictly positive by assumption, this implies the inequality 

e,{x{t)) > {\{tp),R\{tp)) - {A{t),S\{t)) . (9) 

Condition (iii) allows to bound this further to 

eh{x{t)) > (A(tp),i?A(tp)) - {A{t),RA{t)) . (10) 

Due to (i) the time derivative of {A{t),RA{t)) is always non-positive 

dt{A{t),RA{t)) = {A{t), {d R + RC)A{t)) < 0, (11) 

so that the right-hand-side of (10) grows monotonically, which implies that eh{x{t)) > 0 
for t > tp > 0. 
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Whereas x(^p) is (by assumption) positive semi-definite and of rank h, the ultimate 
goal is to show complete positivity for all times f > 0 and to also take into account 
the initial condition e/i(x(0)) = 0. To this end it is necessary to verify that that all 
non-trivial eigenvalues of x(^) become positive in leading order in t. If this is the case, 
and (A(t), RA{t)) is not strictly constant (i.e. it decreases in leading order in t following 
(i)), then (ii) can be replaced by 

(ih) (A(0),(il-5)A(0))=0 


and (i), (ii’), and (hi) verify CP dynamics for all f > 0. 

To prove the latter statement, we first find tp > 0 such that x(fp) ^ 0 ^ind 
eh{x{tp)) > 0, as well as {A{tp), R A{tp)) < (A(0),-R A(0)), and Ai(f) is CP for all 
0 < t < tp. To show complete positivity for t > tp, we rewrite (7) trivially as 

e,{x{t)) = (A(0),5A(0)) - (A(g,RA(g) 

+ {A{tp),RA{tp))-{A{t),SA{t)) . 

Following (ii’) this implies 


e,(x(t)) = (A(0),RA(0)) - (A(g,RA(g) 
+ (A(g,RA(g)-(A(t),5A(t)) , 


(13) 


where the first two terms satisfy (A(0), R A(0)) — {A{tp), R A{tp)) > 0, such that 


enixit)) > {A{tp),RA{tp)) - {A{t),SA{t)) . 


(14) 


This proves strict positivity of eh{xit)) shown above in (9). 


2.2. Geometric interpretation 

The proposed procedure can be intuitively understood in geometrical terms: Just 
like normalized, Hermitian operators are quantum states (i.e. positive semi-definite 
operators) of a two-level system if and only if they are represented by points inside 
the Bloch sphere, also valid extended dynamical maps form a geometrical object 
which we denote by V. Although the dimension of V is substantially larger than 
the dimension of the Bloch sphere, and its shape will typically be more complicated, 
this underlying geometry permits an immediate understanding of our approach. The 
schematic representation of figure 1 depicts all the points within V in blue and yellow. 
In turn, all points outside this area correspond to invalid extended maps with a system 
component that is not CP. The dynamics induced by (4) is indicated by arrows. 

The goal is to find conditions which guarantee that, for a given initial condition 
A{tp) within V (indicated by a black circle in figure 1), the dynamics is such 
that its trajectory will never leave V. As this can typically not be shown 
directly by means of V, one strives for the identification of a region V, such that 

V contains A{t) for all times t > tp 

V lies inside V. 
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Figure 1. A schematic representation of the set V (blue and yellow) of extended maps 
r = 0 \i) with CP Pi is depicted. The arrows show the dynamics according to 

(4). The thick trajectory corresponds to the dynamics that starts at A(0) (black cross) 
and passes through the point A(tp) (black circle). Equipotential lines of the function 
J-(A) = {A,RA) are depicted in green (dotted). The line that contains A{tp) is 
depicted with long dashing and confines the set V (yellow). If V is contained in V (as 
depicted), and iF{A{t)) decreases monotonically, CP dynamics is verihed for t > tp. 

This region V shall be confined by the eqnipotential line of a function -7^(A) = (A, RA) 
(figure 1 depicts some equipotential lines of 7F in green (dotted/dashed)) that is 
determined by R. More precisely, V is given by all extended maps F = XliTj ® |f) 
for which -T(r) < 7F{A(tp)). When the coordinate system is chosen such that the 
asymptotic extended map A{t —)■ oo) lies in its origin, then J^(A(t)) can be considered 
a distance between A{t) and the asymptotic map around which the object V is centred. 

According to (i), the function J^(A(f)) decreases monotonically under the dynamics 
induced by (1), for why J^(A(f)) can never exceed the value 7F[A{tp)) for t > tp. The 
geometric implication of (i) is that A{t) will never leave V. Condition (ii) and (iii), in 
turn, manifest a relation between V and V such that V contains V. This implies that 
and all extended maps inside V are valid. 

The initial extended map A(0) (indicated by a black cross in figure 1) lies on the 
boundary of V (due to x(0) having rank one) for why V and V become tangential in 
the point A(0) when tp becomes infinitesimally small. 

2.3. Simplification of condition (ii’) and (iii) 

The construction of R is typically most challenging when CP dynamics shall be proven 
for f > 0 (i.e. when V and V become tangential in A(0)). A particularly convenient 
case is, however, given when the S is of rank one, which, in geometrical terms, means 
that the valid region V is confined by two parallel planes. The construction of R = R^ 
that proves complete positivity for f > 0 can in this case be simplified as the condition 
(iii) can be partially replaced by a set of equality conditions, which effectively decrease 
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the number of degrees of freedom in R that have to be determined. This is particularly 
helpful for an analytical construction of the latter. 

More specihcally, the inequality R — S > 0 is satished (i.e. (iii) holds true) if R is 
positive semi-dehnite, normalized as in (ii’), and satishes 

^^^(r,i?r)|r=A(o) = 0 (15) 

for a basis of the kernel of S. Equation (15) is a set of linear constraints that 

reduces the number of degrees of freedom in R. 

To understand why this reformulation of (iii) is possible, we emphasize that A(0) 
lies on the boundary of V and can thus not be from the kernel of the operator S. That 
means that A(0) and j qj-q linearly independent and form a complete set (because 
S is of rank one) such that every extended map v can be decomposed as n = aA{0) + bK 
with K being from the span of and a,b gC. Without loss of generality, it can 

be assumed): that a*b G M. With this decomposition, we obtain 

(n, {R -S)v) = a*b{{\{0),RK) + (K, A(0))) + b‘^{K,RK), (16) 

where (ii’) has been employed. From (15), one can, however, deduce 

{A{0), RK) + {K,R A{0)) = lim -[(A(0) + eK, R (A(0) + eK)) - (A(0), A(0))] 

£—>■0 e 

= A(r,iir)|r=A(o) = o, (17) 

which, together with R > 0, implies non-negativity of (16) and thus positive semi- 
dehniteness of i? — S'. 

With the framework derived hitherto, we are now in the position to analyse 
exemplary equations of motion that permit to construct R by purely analytic means. 
Subsequently, in section 3 we will expand the framework of constructing R towards 
numerical techniques that will permit to treat more complex cases as exemplihed in 
section 3.3 and 4.3 

2.4- Examples 

2.4- 1- The damped Jaynes-Cummings model An instructive example of a HEOM is 
given by the damped resonant Jaynes-Cummings model, which describes the dynamics 
of a two-level system inside a leaky cavity. The electro-magnetic cavity held is 
characterized by a Lorentzian spectral density function that is centred around the 
transition frequency of the two-level system. The spectral width of the cavity held 
is denoted by (, whereas 7 labels the coupling strength of the two-level system and the 
cavity held. After tracing out the electromagnetic held, this model permits an exact 
description of the (generally non-Markovian) dynamics in terms of a HEOM, which reads 

f If that is not the case, then one can consider b' = 6-|-5(a*&)/9(a) and K' = (1 —3(a*6)/[6'3(a)])/ir G 
ker(ii) such that V = aA(0) + b'K' and ^{a*b') = 0 
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Qi{t) 

= C22(f) 

(18a) 

22 (t) 

= 7^^-2i(f) +C X] criP2(f)cr| +C23(f) 

(186) 


i=x,y,z 


23 (f) 

= 1 X] 22 (f) -2C alg3{t)a. 

(18c) 


i=x,y 


with.V_ Qi = —{cr_|_(j_, ^)i}/2 and cr± = {ax^iay)/2. The system state is denoted 

by Qiit), whereas the two anxiliary operators Q 2 {t) and Qzit) encode information abont 
the environment and system-environment correlation. The introdnction of a general 
framework by means of which this HEOM has be obtained is given in Appendix A.l and 
remarks that explicitly concern the constrnction of (18a) - (18c) are fonnd in Appendix 
A.l.3. 


To tnrn this HEOM eqnation into an eqnation of motion for the dynamical maps 
Ai(t) to A 3 (f), one writes (18 a)- (18 c) in the Bloch basis. The Bloch representation of 
the extended state is then given by a 12 -dimensional vector, and the 12 x 12 -dimensional 
Bloch representation for the generator C, is also the generator for the dynamics of the 
12 X 4-dimensional Bloch representation of the extended dynamical map. The latter is 
initialized by Ai(0) = I and Aj(0) = O for i > 2, where I and O are the identity and the 
nnll map, respectively. This corresponds to the initial condition P2(0) = ^ 3 ( 0 ) = O. It 
has been algebraically verihed that the coefficient matrix x{t)^ which characterizes the 
system map Ai (see (5) with /ij = cTj), is always of the form 


x{t) = 


1 

4 


(Ai(f) + 1)2 0 0 

0 1 — A^(f) i — iXfit) 

0 ^Xl{t) — i 1 — Xl{t) 
Xl{t) -1 0 0 


m -1 
0 
0 

(Ai(f)-1)2 


(19) 


where the only time-dependent degree of freedom Ai(t) is a coordinate of the system map 
Ai(f). Thns, the condition x(f) > 0, i.e. complete positivity of the system map Ai(f), 
does only depend on Xi{t) and all other degrees of freedom of Ai{t) can safely be ignored. 
The dynamics of Ai(f) is obtained from the differential eqnation dtA{t) = C A{t) and 
reads dtX(t) = I X(t), where X(t) = [Ai(t), A 2 (t)] is a two-dimensional real-valned vector 
and 


I 


0 C 

-7 

2 


( 20 ) 


The dynamical variable Xi{t) is a coordinate of Ai{t) and the initial condition A(0) 
corresponds to A(0) = [1,0]. Strict positivity of all non-trivial eigenvalnes of x(f) (two 
of them vanish constantly) in the first non-vanishing time-step is given when xC > 0 . 
When this prodnct is positive, we can replace the normalization condition (ii) by (ii’). 
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Figure 2. Figure (a) depicts a representation of the set V (blue/yellow, solid) of 
extended dynamical maps with a system component that is CP. The green (dashed) 
equipotential line of the monotone E contains the initial condition A(0) (black cross) 
and bounds the set V (represented in yellow). The dynamics induced by (18a)-(18c), 
is depicted by vectors (grey/black) and the trajectory for the initial condition A(0) is 
indicated in black (thick line). The set V lies completely inside V and as E decreases 
monotonically, this proves complete positivity. Figure (b) shows eh{x) as a function 
of time and gives an idea of what values the smallest non-trivial eigenvalues of x(<) 
can take. Blue (solid) corresponds to 7 = 10^ (same parameters as in (a)), whereas 
red (dashed) represents 7 = 10®C- Both cases are covered by the conditions derived in 
section 2.4.1. 


Due to the structure of xif) (19), its highest and the second-highest elementary 
symmetric polynomials e^i^xit)) ^-nd e^^xit)) cire constantly equal to zero. The complete 
positivity of the system map Ai(f) thus depends on the highest non-vanishing polynomial 
^h{x) = 62 (x) = — tr(x^))/2 = (1 — Af)/4 which is positive if and only if 

|Ai(f)| < 1. This condition is expressed as eh{x{t)) = A'*'(0) s A(0) — A'^(t) sX{t) > 0 with 
a matrix s = diag[l, 0]/4 that characterizes the valid region V. The situation is depicted 
in hgure 2(a). 

To prove strict positivity of e/i(x(f)) for all times t > 0, we strive for a matrix r 
(the representation of R in the space of X{t)) that satishes r + r I <0 (condition (i)), 
A'^(O) (r —s) A(0) = 0 (condition (ii’)), and r — s > 0 (condition (hi)). Condition (ii’) and 
(iii) hold true, when r is parametrized by r = diag(l, r 22)/4 with r 22 > 0. Condition (i), 
in turn, is satished if g = /I r -|- r / < 0 , which is the case if and only if trace and 
determinant of —q are non-negative. The latter is given as det(—g) = —( 7^22 — 2(C)^/4 
for why we have to choose r 22 = 2^/7 and with tr(—g) = 2 r 22 C = H fh® conditions 
g < 0 and r 22 > 0 are equivalent to 7 > 0 and C > 0. Together with 7 (C > 0 this proves 
that positivity of both parameters 7 and A is sufficient for eh{x(f)) > 0 and implies 
CP dynamics. Even though eigenvalues can become very small (see hgure 2(b) for an 
example), the conditions clearly assert that no eigenvalue of x(f) vanishes. 
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2 . 4 . 2 . Reviving coherences: a HEOM with two levels Another model for a non¬ 
monotonic decay of phase coherence can be formnlated in terms of the hierarchy eqnation 
(see Appendix A. 1.1 for an explicit derivation) 

Pi(f) = +ujg2{t) (21a) 

g2(t) = aojD^Qiit) +^2(TzQ2(t)az {21b) 


with the Panli operator cr^ and the dephasing Lindbladian D^g = Czgc^z — g. For 
a = 0 and the initial condition ^ 2 ( 0 ) = 0, this hierarchy eqnation describes pnrely 
Markovian dynamics with exponentially decaying coherence. For non-vanishing values 
of a, however, the dynamics deviates from a mere exponential decay and can give rise 
to genuine non-Markovian revivals. The initial condition ^ 2 ( 0 ) governs the infinitesimal 
behaviour of gi{t) at time t = 0. A generic choice is ^ 2 ( 0 ) = 0 although other choices 
are possible as we will demonstrate later. 

In the Bloch basis, the extended state q can be expressed in terms of an eight¬ 
dimensional vector and the according representation of the generator jC reads 


0 

0 

0 

0 

U) 

0 

0 

0 

0 

-7i 

0 

0 

0 

OJ 

0 

0 

0 

0 

-71 

0 

0 

0 

UJ 

0 

0 

0 

0 

0 

0 

0 

0 

UJ 

0 

0 

0 

0 

72 

0 

0 

0 

0 

—2au 

0 

0 

0 

-72 

0 

0 

0 

0 

— 20 : 0 ; 

0 

0 

0 

-72 

0 

0 

0 

0 

0 

0 

0 

0 

72 


With the initial conditions Ai(0) = I, and A 2 ( 0 ) = O and the structure of C one 
can show by purely algebraic means that the Bloch representations for the extended 
dynamical map \{t) reads diag[l, Ai(f), Ai(f), 1] (g) |1) -|- diag[0, A 2 (t), A 2 (t), 0] ( 8 ) | 2 ) 
and depends only on the two scalar parameters Xi{t) and X2{t). With this form 
of the dynamical map Ai(f), the corresponding matrix x{t) is given by x(f) = 
diag[l -|- Ai(f), 0,0,1 — Ai(f)]/2. In particular, the maximal rank of x(f) is h = 2 
and complete positivity of Ai(f) is given for |Ai(f)| < 1, or, equivalently, e/i(x(t)) = 
Ai(0)sA(0) — X^{t)sX{t) > 0 with s = diag[l,0]. Strict positivity of all non-trivial 
eigenvalues of x(f) for infinitesimal times is given when either 71 > 0 , or 71 = 0 and 
u},a > 0 hold true. 

The procedure of proving complete positivity is now analogous to the previous case: 
As Ai(f) and X2{t) are the only relevant dynamical variables, the equation of motion for 
the extended dynamical map can be reduced to dtX{t) = I X{t) with X{t) = [Ai(f), A 2 (t)], 
the initial condition A(0) = [1,0] and the generator 


I 


-71 u 
—2au —72 


(23) 
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Figure 3. Time evolution of the coherence tT:{(JxQi) for an initial state 0 i(O) = |+)(+| 
with 1+) = (|0) + |l))/-\/2, evolving according to (21a), (216), and its extension (25). 
The dynamics of (21a) and (216) with initial conditions 02 ( 0 ) = 0 and 0 i(O) = 0 is 
depicted in blue (solid) and red (long dashing), respectively (71 = 72 = l/2a;, a = 4). 
The orange (short dashing) line represents the dynamics induced by introducing the 
additional level in (25) with ^ 2 ( 0 ) = ^ 3 ( 0 ) = 0 (/3 = 8, 7 = l/2a;). 


As already seen for the Jaynes-Cummings model, the operator R is characterized 
by a 2 X 2-dimensional matrix r that is parametrized by r = diag[l, r 22]/4 with r 22 > 0 
in order to satisfy the conditions (ii’) and (iii). The monotonicity condition (i) holds 
true if /I r -h r Z < 0 , and a good choice for r 22 (the choice is not unique) is given by 
f 22 = {auj^ + 7 i 72 )/( 2 a^ci;^). In this case, monotonicity and r 22 > 0 are equivalent 
to 7 i > 0 , 72 > 0 and 2 a + 7172 > 0 such that these three conditions together with 
^ 2 ( 0 ) = O and strict positivity of x(f) for short times are sufficient for the hierarchy in 
(21a) and (216) to induce CP dynamics. An example of a process that satishes these 
conditions is given in hgure 3. 

It is well established that the initially linear decay of coherence, that is depicted 
in the inset of hgure 3, is an artefact of the Markov approximation. In the framework 
of hierarchy equations, it is straight-forward to remove this artefact by choosing the 
initial condition ^ 2 ( 0 ) = —£iiPi(0), which corresponds to A(0) = [ 1 , 71 /a;] and yields 
Pi(0) = 0. Figure 3 shows the difference of the time evolution in the hrst time-steps. 

The examination of complete positivity requires a modihcation of the matrix r in 
order to still comply with (ii’) due to the change of the initial condition A(0). To this 
end, we parametrize r by means of the Cholesky decomposition such that it is positive 
semi-dehnite and A(0) (r — s)A(O) = 0, which (due to the symmetry of r) reduces the 
number of degrees of freedom in r to two. The condition r — s > 0 is then incorporated 
into the parametrization of r as shown in section 2.3, which eliminates another degree 
of freedom. With this parametrization, the two conditions (ii’) and (iii) are satished 
and all but one single degree of freedom of r have been determined. This last degree of 
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freedom can for example be chosen snch that det(/f r + r /) = 0. In this case one obtains 


^ 7? + ^ -W7i 

4:7] —cvyi 


(24) 


with Tj = ^auP' + 7 i 72 . The monotonicity condition (i) and r > 0 are then eqnivalent 
to 7 i + 72 > 0 and 2 aa;^ + 7172 > 0 , which (together with the condition that all non¬ 
trivial eigenvalnes of x{t) become strictly positive for short times t) are thns snfficient 
conditions for CP dynamics. 


2 . 4 . 3 . Reviving coherences: a HEOM with three levels The HEOM in (21a) and (216) 
can be extended by adding a new anxiliary operator ojgsit) to the right-hand side of 
(216) and angmenting the eqnation of motion by (see Appendix A. 1.2 for a rigorons 
derivation of this eqnation) 

hit) = I3uj(^z 02 it) (tI + 73 Qsit) 0 -^ • (25) 


Snch HEOM gives rise to a broader variety of dynamical processes and rekindles the 
qnestion of complete positivity. In favonr of a concise discnssion, we rednce the nnmber 
of free parameters by setting 71 = 72 = 73 = 7 . 

The extended state git) as well as the generator C are again expressed in the 
Bloch basis. In similarity to (22), this permits to express the HEOM dtQ{t) = JC Q{t) 
in terms of a vector eqnation and the transition from states to dynamical maps is 
made by replacing git) by the extended dynamical map A(t), which is initialized by 
Ai(0) = I and Aj(0) = O for i = 2,3. In the Bloch representation, all dynamical 
maps Aj( 6 ) preserve the strnctnre diag[hii, Ai( 6 ), Ai( 6 ), ha] snch that the dynamics of 
A(6) is completely characterized by the three-dimensional real-valned vector A(6), which 
is initialized by A(0) = [1,0, 0] and evolves according to c?tA(6) = I A(6) with 


I 


—7 u 0 
— 2 au —7 u 
0 —/So; —7 


(26) 


In fact, the qnestion of complete positivity only depends on the two ratios a = auh 
and /3 = as one can see in terms of the scaled variables Aj( 6 ) = Aj( 6 )a;Y 7 b 

whose eqnations of motion is described by I that is obtained from (26) by mnltiplying 
all elements kj with (a;/ 7 )*“T Since the strnctnre of 7 ( 6 ) does not change as compared 
to the HEOM with two levels, its highest non-trivial elementary symmetric polynomial 
= 62 is still correctly expressed by Chixit)) = Af(0)sA(0) — Af(6)sA(6) with 
s = diag[l, 0, 0]/4. Again, a Cholesky decomposition of r (snch that r > 0) and the 
condition Af(0)(r — s)A(O) = 0 permit to express r —s > 0 in terms of two linear eqnality 
constraints (see section 2.3), which give rise to the parametrization 


r = 


1 

4 


1 0 0 

0 a ca 

0 ca ac^ + b 


(27) 
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Figure 4. Parameter set for which the system dynamics induced by the HEOM in 
(21a), {21b) and (25) is found to be CP. Green (short-dashed) corresponds to purely 
analytic verification of complete positivity, and orange (long-dashed) is based on a 
numerical optimization of the operator R. White corresponds to violation of complete 
positivity identified through explicit solution of the equations of motion. The blue 
region (solid), thus seems to correspond to CP dynamics, but can not be identified as 
such with the present framework. 


with a,b > 0. Considering the more interesting case of d, /3 7 ^ 0, a very convenient 
choice for the three free parameters is given by a = 1/|2 q!|, b = l/|2d/3| and c = 0 . 
This matrix satisfies the monotonicity condition Vr -|- rZ < 0 if and only if one of the 
following two conditions is satisfied: 

— 1 / 2 <q ;<0 and —2a — l<j3 (28a) 

d > 0 and —1 </3 . {28b) 

Conseqnently, any of the two conditions (28a) or {28b) ensnres complete positivity when 
strict positivity of all non-trivial eigenvalnes of x{t) is given for short times. 

Althongh this is only one instance of how the matrix r can be chosen and other 
choices (e.g. with c 7 ^ 0 ) are possible, it can prove complete positivity for a rather large 
parameter regime as shown in fignre 4. An individnal nnmeric optimization over r for 
each nnmerical instance of the HEOM can even slightly enlarge that parameter set. 

3. Numeric approaches to more complex problems 
3.1. Transformation of the extended dynamical map 

A central element of the procednre discnssed in section 2 is to express the highest 
non-trivial elementary symmetric polynomial eh{x{t)) in terms of a qnadratic fnnction 
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in A(t) as shown in (7). For the examples discussed in section 2.4, such a representation 
has been obtained by anticipating the structure of the matrix x(t) that characterizes 
potential solution to the initial value problem by purely algebraic means. Because the 
initial condition Ai(0) = I and Ai(0) = O for i > 2 is fixed, it is for example often 
possible to prove that particular subspaces are never occupied. To this end, one shows 
that the dimension of the space that is spanned by C ^ A(0) with /c = 0,..., d and d being 
the number of degrees of freedom in A{t) is strictly smaller than d. When that is the 
case, then this indicates the existence of preserved quantities, which permit predictions 
about the form of x{t) ^md can thus help to obtain (7) (and to reduce the dimensionality 
of A(f)). 

To make the approach more universal and to permit eh{x{t)) to be of a degree 
m > 2 in x{t)j can apply a transformation which renders geometries that are always 
accessible to our technique. To this end, we define a function 

G[x{t)] =c^ - (ehixit)) - cf = 2ceh{x{t)) - el{x{t)) 

< 2ce/*(x(t)) , 

with a constant c > 0 such that the strict positivity of is sufficient for 

Chixit)) > 0. Although G[x(t)] is of a higher degree in A(f), considering the function 
rather than Chixit)) has the big advantage that the prior can always be expressed 
in a form that is very similar to (7). This is done by defining a new vector S(f) whose 
first component reads Si(f) = Chixif)) Independently of the components Sj(f) with 
i > 2, this vector permits to express as 

G|xM| = 5t(0)55(0) - 5t(t)gg(t) , (30) 

with a matrix S that is of the simple form S = diag[l,0,... ,0]. It is essential that 
E(t) follows a linear differential equation dfEit) = CE{t) and as the elements of E(t) 
are functions of x(f) (^-nd thus A(t)), this differential equation needs to be derived 
from the HEOM. Already the first component Si(f), however, is a non-linear function 
in A{t) and its time-derivative does typically also involve terms that are non-linear in 
A{t). For H(f) to anyway evolve according to a linear equation, one can now consider 
any non-linear contribution to be a new element of H(f) and keep increasing the size of 
S(f) until one has recovered a linear equation of motion cltS(f) = §. The matrix 

C is constructed accordingly and eventually yields the desired equation of motion for 

m- 

In order to prove strict positivity of G[x{t)], one follows the same procedure that 
has been introduced in section 2.1 and seeks a matrix R > 0 that satisfies (i) to (iii) 
(where C and A{t) must be replaced by £ and E{t) and the scalar product becomes the 

§ A simple example would be the scalar equation of motion x{t) = ax(t) and the non-linear coordinate 
transformation y{t) = x^{t)+x{t). The coordinate y{t) satisfies y{t) = 2ay{t) — ax{t); replacing x{t) in 
terms of y{t) would result in a non-linear equation of motion, but treating x{t) as independent variable, 
permits to maintain linear equation of motion at the expense of increasing the number of independent 
variables. 
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standard scalar product in the Euclidean space). To this end, it is beneficial to choose 
c = linii^oo ^hixif)) such that the second term in (30) vanishes asymptotically. If such 
a matrix R can be found and G[x(t)] increases in first order, then G[x{t)] > 0 hold true 
for t > 0 and the system dynamics is completely positive. 

3.2. Semi-definite programming 

The transformation introduced in section 3.1, or a consideration of more complex 
quantum systems and HEOMs with several levels, can quickly increase the 
dimensionality of the HEOM to a level at which the operator R can not be constructed 
analytically. In such cases, however, the convexity of the problem permits the 
application of a highly efficient numerical approach. To this end, let us parametrize R 
such that the normalization condition (ii) (or (ii’)) is satisfied and define Q = R+RC 
and B = (nl — Q) © {R — S). The minimization problem 

Vm = min (n I .B > 0) (31) 

v,R 

is a semi-definite program [29], which can be solved reliably and efficiently. The 
constraint B > 0 implies R — S > t) (i.e. condition (iii)) and Q < nl. If the minimum 
Vm is non-positive, one has found R such that (i) holds true. The condition Um < 0 thus 
verifies complete positivity. 

3.3. Example: non-Markovian dynamics due to a finite temperature bath 

To gain a better intuition of how the transformation introduced in section 3.1 affects 
the geometric properties of the problem, we consider a two-level system that is coupled 
to a bath of finite temperature. In the Markovian case, the situation is described by a 
master equation in Lindblad form = B+gi B-Qi where 

B±gi = 7± /2) , (32) 

with rates 7 ±. In order to also permit non-Markovian dynamics, we introduce the 
traceless auxiliary operator gi{t) and extend the Lindblad equation to the HEOM (see 
Appendix A .2 for a motivation of this equation) 

dtgfit) = {B+ + B_)gi{t) +Ujg 2 {t) (33a) 

dtg 2 {t) = —iB+ + B.) gfit) (P. + Vy) g^it) (336) 

7 + 2 

with 7p = (7+ + 7-)/2, 7m = (7+ — 7-)/2, and ^ and u being scalar parameters that 
characterize the additional levels in HEOM (the special cases ca = 0 and = 0 give 
rise to Markovian dynamics). With g = pi © |1) + p 2 © |2) and the operator jC that 
is derived from the HEOM (33 a)- (33 6), the latter can be expressed as digit) = C g(t). 
The equation of motion for the dynamical maps is obtained by replacing g{t) with 
A = Ai © |1) + A 2 © |2) that is initialized by Ai(0) = I and A 2 ( 0 ) = O. 
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The structure of C and the initial condition give rise to a set of identities and 
restrictions which hnally permit to fully represent A(f) by means of a four-dimensional 
vector A(t), whose initial condition is A(0) = [— 7 ^/ 7 ^, 0,1,0], and which evolves 
according to dtX(t) = I X(t) with 


27 p —u 0 0 

2e 27 p 0 0 

0 0 7p —oj 

0 0 ^ 7p 


(34) 


The dynamical map Ai(f) in the Panli basis is given as in (5) (with /Xj = ai) and 
dne to the strnctnral properties of C and A(0), the coefficient matrix is always of the 
form 



27m, A 3 (t)+7m -7p1 (b 

0 

0 

e(t) 1 


47m 

47p 


0 

m 


0 

x{t) = 

0 

47m 

i6{t) 

47p 

m 

0 


47p 

47m 


m 

0 

0 

-27mA3(b-|-7m-7pAl(t) 


47p 

47m 


(35) 


where 9(t) = 7 ^ + 7 pAi(f). In contrast to the previons examples, this strnctnre of x(f) 
does not permit conclusions abont constantly vanishing eigenvalnes snch that the highest 
possible elementary symmetric polynomial Ch = det mnst be considered. Eqnation (35) 
does, however, give rise to the factorisation Chixif)) = det(x(t)) = aPf{t) P 2 it) with a 
proportionality factor a = 7 + 7 _/[ 4 ( 7 ^ — 7 +)]^ and 


-Pi(^) = IpXiit) + 7m (36 a) 

P2{t) = Cl — C 2 X'^it) C 3 Ai(t) -\- C 4 A^(t) . (366) 

Here, the constants are given by Ci = 7 _ 7 + 7 ^, C 2 = 7 ^ — 7 +, C 3 = 7 ! — 7 ^ and 
C 4 = 7-7+7p- 

Considering the case 7+, 7- > 0, the proportionality factor a is positive snch that 
it is sufficient to show strict positivity of Pf{t) and P 2 {t) separately in order to prove 
the positivity of eh{x{t))- After verifying strict positivity of x(f) for short times, this 
ensnres CP dynamics for f > 0. 

The factor Pf{t) is always non-negative bnt to show strict positivity for f > 0, we 
reformnlate the problem snch that Al(0) A(0) — Al(t) X{t) > 0 with a symmetric 

matrix is sufficient for Pf{X{t)) > 0 and strive for the constrnction of a matrix 
r^^'> that satishes r^^^ + r^^'> I < 0, Al(0) A(0) = 0, and > 0 (see 

(i), (ii’), and (iii)). The procednre is very similar to the previons examples for why we 
abstain from a detailed discnssion and anticipate that Pi{t) never retnrns to zero when 
— 27 p < oj^ holds trne. 

Showing the second condition P 2 {X{t)) = Al(0) A(0) — Al(t) A(t) > 0 is 

more involved, which can be nnderstood best in geometric terms. A cross-section (for 
A 2 = A 4 = 0) of the set V 2 of vectors A that satisfy P 2 {X) > 0 is depicted in hgnre 3.3(a) 
and the black cross corresponds to the initial condition A( 0 ). As one can see, it is not 
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(a) 

0.05 

Ai 0 

-0.05 

- 0.1 


Figure 5. The sets inside which P 2 > 0 is satisfied are depicted in blue before (a) 
and after (b) the non-linear transformation from X{t) to S(t). The black lines show 
the trajectories of X(t) and S(t) for the initial conditions that correspond to Ai(0) = I 
(black crosses). For the geometry depicted in (a), it is impossible to find a matrix 
such that the equipotential line of points ^ with ^ ^ = const contains A(0) and 

lies completely inside V 2 (the set with P 2 >0). After the non-linear transformation a 
more benign geometry has been obtained such that an according matrix can be found. 



possible to place an ellipsoidal equipotential line A = const (that would confine the 

set V 2 ) inside V 2 such that it contains the initial condition A(0). In algebraic terms, the 
two conditions (ii’) and (iii) can not be satisfied simultaneously for geometric reasons. 

We can, however, transform the geometry of V 2 in order to mitigate this problem. 
To this end, we define the first component of the new vector S(f) to be given by 
Si(t) = c — P 2 {t) with c = Iim 4 ^oo(-P 2 (^)) > 0 such that G = Sl(0)^S(0) — 
E'l{t)SE{t) > 0 implies P 2 {t) > 0, when S = diag[l, 0,..., 0]. The other elements 
= /j(Afc(f), Afc(f) Xiit)) with j > 1 are polynomials in Afc(f) of degree one or two and 
are iteratively chosen such that E(t) satisfies the differential equation dtE{t) = 
with an according generator C. 

Figure 3.3(b) depicts the convex set of vectors S for which G > 0 holds true. 
Because of its benign geometry, we can find a matrix R such that the equipotential line 
Ef RE = const contains the initial condition S(0) and lies completely within V 2 , i.e. 
Sl(0) {R — S') S(0) = 0 and R — S > 0 hold true. When also the monotonicity condition 
i? -b i?G < 0 is satisfied, then this proves P 2 (t) > 0 for all times. 

Due to the rather high dimensionality of E{t), such a function can not be 
found analytically but by means of a semi-definite program a very efficient numeric 
construction has been carried out. 

4. Quantifying the violation of complete positivity 

In many cases, an exact description of the dynamics of a physical system is only 
obtained for infinitely deep hierarchy equations. A truncation of the HEOM is then often 
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employed in order to obtain approximate system dynamics. Such a truncation, however, 
can give rise to a slight violation of complete positivity of the induced dynamical system 
map Ai(f). The latter typically occurs for short times, at which the coefficient matrix 
x{t) has very small eigenvalues, and vanishes for all times greater than a critical time 
tp. During this initial time interval 0 < t < tp complete positivity is violated and even 
if the violation is numerically negligible, the proposed framework will not be able to 
verify CP dynamics. 

In such cases, it is desirable to estimate the degree of violation of complete positivity. 
This degree is for instance an indicator of whether the approximation is good enough 
as strong violation of complete positivity suggests unphysical dynamics and thus an 
insufficient approximation. We will propose two different ways of quantifying this 
violation and demonstrate their applicability by means of an important example. 

4 . 1 . Complete positivity after an initial violation time 

At the initial time f = 0, when the system map Ai(0) is equal to the identity, the 
coefficient matrix x(0) has just one non-vanishing eigenvalue. A truncation (and thus 
approximation) of the HEOM therefore often causes an eigenvalue of x{t) to become 
negative for short times. The identification of a critical time tp at which the dynamics 
becomes CP again helps to estimate whether the violation is indeed just a short-time 
phenomenon, or a structural problem that eventually impedes a physical interpretation 
of the time-evolution. A good guess for tp can be obtained from a numeric propagation 
of xit) foi' short times and once a candidate for tp (and the according extended map 
A(fp)) has been identified, the proposed method in section 2.1 for proving complete 
positivity can be employed to show that the dynamics is CP for all t > tp. 

The situation is depicted in figure 6(a), where it is schematically visualized how the 
extended dynamical map leaves the valid region V for short times in order to re-enter 
at time fp > 0. 

4 . 2 . Lower bound on the eigenvalues of x 

An alternative way of quantifying the violation of complete positivity of the system map 
Ai(f) is based on estimating a lower bound —5mm on the eigenvalues of the matrix x(f). 
Rather than requiring the positivity of the highest non-trivial elementary symmetric 
polynomial e/i(x(t)), we strive to prove Chixif) + ^mirJ-) > 0, which permits negative 
eigenvalues of x{t) as long as they remain greater than the new parameter —Smin- 

Algebraically this means that one needs to find an extended map El = C) |i) 

that is such that (7) is modified to 

eh(x(t) + 5m*nI) = (f^,5Pl)-(A(f),5A(f)) . (37) 

Furthermore, the normalization condition (ii) is then replaced by 

(ii”) {n,{R-s)n) = 0 

and the set of conditions on R is extended by 
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Figure 6. The figure shows schematically how a time-evolution violates complete 
positivity for short times before the dynamics becomes CP again. The induced 
dynamics that is initialized by A(0) (black cross) is represented by the black trajectory. 
The valid region V of extended maps that induce CP dynamics is depicted in blue and 
A{t) leaves this region for short times. In (a), the initial condition is thus replaced 
by A{tp) (red cross). The time interval 0 < t < tp is disregarded and the proposed 
procedure is employed to determine the set V, that contains A(t) for t > tp, such 
that V lies inside V and complete positivity is shown for t > tp. In (b), the condition 
xit) > 0, that defines V, is replaced by x(t) > —6min^ which enlarges the set V of 
states that are considered valid. The tangential point Q oiV and V becomes object to 
an optimization, which increases the chance to find V such that V CV. In that case, 
—Smin is a lower bound for the eigenvalues of x(t). 


(iv) {n,Rn) > (A(o),i^A(o)). 

Whenever R is found such that (i), (ii”), (iii), and (iv) are satished, this proves strict 
positivity of Chixif + ^min)) for all times t > 0. 

To prove this statement, (ii”) is applied to (37), which yields 

eh{x{t) + = {^,Rn) - {X{t),S X{t)) . (38) 

By means of (iii), this can be further bound to 

ehixit) + ^mirJ) > (Xt,RQ,) - {X{t),R\{t)) , (39) 

which, according to (iv), is strictly positive at the time t = 0. For all times t > 0, the 
right-hand side of (39) increases monotonically due to (i), such that Chixif) + ^mirJ) > 0 
and all eigenvalues of xit) are greater than —5min for all times. 

In geometric terms, this modihcation enlarges the region V that is considered valid 
as depicted in hgure 6(b). The points F with J^(r) < -F(r2) (where is dehned as in 
section 2.2) dehne the region V that will contain the solution A(t) for all times and El is 
the tangential point of V and the (enlarged) set V. It is important to point out that f2 
is not uniquely dehned and becomes itself object to an optimization. This often permits 
to hnd an operator R that satishes (i), (h”), (iii), and (iv) and proves —hmin to be a 
lower bound on the eigenvalues of xit)- 
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4 . 3 . Example: The spin-Boson model 

One of the most prominent examples of a microscopically derived HEOM describes the 
dynamics of a two-level system that is conpled to a bath of harmonic oscillators. For 
a high temperature case (/dy <C 1) with weak system-bath coupling (A^/y^ cn), the 
HEOM is expressed as [10, 12, 13] 

Qi = -iA[cr„,P2] (40 a) 

A/3y oj 

Q2 = -iA[(T^,Qi] - -i^Wz,g2]-702 ■ (406) 

Here, gi{t) represents the state of a system with a resonance frequency u, whereas [•, ■] 
and {■, ■} denote the commutator and the anti-commutator, respectively. As discussed 
in [10, 12, 13], A characterizes the system-environment coupling, f3 corresponds to the 
inverse bath temperature, and y is a damping constant. 

Again (40a)/(406) can be written as dtQ{t) = C Q{t), where ^ 0 |1) -|- p 2 ® |2) 

and the generator C has been deduced form the HEOM. In order to obtain the 
equation of motion for the extended dynamical map A(f), one replaces Q{t) by A(f) = 
Ai (g) |1) -|-A 2 0 |2). Expanding the eigenvalues Ciit) of the matrix xif) that characterizes 
the dynamical map Ai(t) for the system state in time, one hnds that one of these 
eigenvalues is given by ei(x(f)) = (/3^y^A^a;^) + 0{t^). For non-vanishing values 

of ft, y, A and u, this eigenvalue becomes negative for short times, which results in a 
violation of complete positivity of Ai(f). Sufficient conditions for complete positivity 
can hence not be found for this important example. It is, however, possible to estimate 
the degree of violation of complete positivity: A numerical integration yields that the 
violation is typically not very strong and limited to rather short times, both of which 
we prove in the following. 

4-3.1. Complete positivity for t > tp By propagating the extended dynamical map 
over a short period of time, it is possible to determine a time tp at which the coefficient 
matrix xitp) is positive dehnite. Based on the procedure laid out in section 4.1, we want 
to prove that complete positivity is given for all subsequent times. 

Due to the structure of C and A(0), it is possible to represent the extended 
dynamical map A{t) by six degrees of freedom that form the six-dimensional real-valued 
vector X(t), which is initialized by A(0) = [1, 0,1,0, 0,1]. The latter evolves according 


to the differential equation dtX{t) 

= l 

X{t) 

with 




■ 0 

u 

0 

0 

0 

0 



0 

0 

—u 

0 

0 

0 


1 = 

7 

0 

4- OJ 

U) 

0 

-7 

0 

0 

-7 

0 

—u 

0 

0 

(41) 
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0 
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The simplification of A{t) yields a set of identities for the dynamical map Ai(t), which 
gives rise to in a certain structure of x(f) given by 


x(*) = 


Xii(^) 

0 

0 

Xu{t) 

0 

X22{t) 

X23(^) 

0 

0 

A23(^) 

Xssit) 

0 

xUt) 

0 

0 

X44it) _ 

ats are 

linear functions 

in X{t) 


(42) 


permit conclusions concerning the existence of constantly vanishing eigenvalues such 
that the determinant is considered to be the highest non-trivial elementary symmetric 
polynomial. It is evaluated to Chixif)) = det(x(t)) cx Pi{t)P 2 {t) with a positive 
proportionality factor and 


Fi = - 1)A4 - (Ag - 2)A6 - 4) - 4(Ai - X^f + b (43a) 

P2 = -4 {AXl + (Ai + As)' - (As + 1)2) - (2/3AA4 + - 1))' , (436) 

where b = + 4. For the considered initial condition A(0), both factors vanish 

initially and by algebraic means it can be shown that Pi{t) = P 2 {t) holds indeed true 
for all times f > 0. It is thus sufficient to keep track of Pi{t) only, and to show that 
it never returns to zero, in order to prove e/i(x(t)) > 0 and thus complete positivity of 
Ai(f) for t > tp. 

For geometrical reasons that are similar to those discussed in section 3.3, the 
transformation that has been introduced in section 3.1 is required in order to render a 
situation in which our procedure is applicable. To this end, a vector E{t) is constructed 
based on X{t) that is such that G = Sl(0) AS(0) — Sl(t) S'E{t) > 0 (with the matrix 
S = diag[l, 0,..., 0]) implies Pi(t) > 0 for t > tp. The strict positivity of G(t) is proven 
when a matrix R has been found that satisfies C"' R + RC < 0, Sl(fp) {R — S) S(fp) = 0, 
and R — S > 0, which can be carried out by means of a semi-definite program. 

As an example, we consider the case 7 = 3a;, A = 2a; and ft = 8/10a;“^ depicted 
in figure 7. This case features a clear signature of non-Markovianity, which is measured 
[30] to A/'(Ai(t —)■ CX))) = 0.34. Complete positivity is proven after the critical time 
tp ~ 0.6/a;. It is noteworthy that the non-Markovian revivals arise after the critical 
time tp. 


4 . 3 . 2 . Lower bound on the eigenvalues of xif) The second possibility to quantify the 
violation of complete positivity is to estimate a lower bound —Smin on the eigenvalues 
of the matrix xif) as discussed in section 4.2. To this end, we make a guess for 6min > 0 
and prove euixif) + bmirJ) = det(x(f) + 5mirX) > 0 for t > 0. The time-evolution of Aft) 
is not affected by this modification of Ch for why the generator I given in (41) and the 
structure of the matrix xft) in (42) still apply. 

The modified polynomial Ch can again be factorized to eh{xft) + oc P[ft)P 2 ft) 

with a positive proportionality constant and 
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Figure 7. A measure of non-Markovianity [30] A/'(Ai(f)) is depicted as a function of 
time for two different numerical instances of the HEOM in (40 a) and (406). The blue 
dots represent the non-Markovianity for a HEOM determined by 7 = 3a;, A = 2a> 
and /3 = 8/10a;“^, which converges to a non-Markovianity of —>■ 00 )) — 0.34. 

Our framework proves the dynamics to be CP for all times greater then tp ~ 0.6/a;. 
The red triangles relate to the case of higher temperatures (7 = lu), A = 2/10 a; 
and /3 = 2/10 0 ;“^), for which the total non-Markovianity is measured to N{Ki{t 
00 )) = 0.08. As we show, none of the eigenvalues of the matrix x(t) ever falls below 
—Smin = —10“^ such that CP-violation is moderate. 


~ + 325mm (1 + ~ -^ 6 (^)) (44a) 

^2^^) — P'lit) + 325mm(l + 25mm + •^ 6 (^)) ) (445) 

where Pi, 2 it) are given in (43a) and (435). In contrast to the previous case, the two 
factors P[ 2 {'t) differ from each other for why they must be examined individually. 

The procedure to prove P[ 2 (f) > 0 is similar to the previous case in section 4.3.1. 
Again, transformations are applied that give rise to new vector representations 
E^^^t) and S(^)(f) for the extended dynamical map and the according generators 
(with i = 1,2) which govern the dynamics. Matrices and functions Gi{t) = 

2 h)(t) are dehned such that Gi{t) > 0 implies T)'(t) > 0 
for all times. The identification of matrices that satisfy [£(*)]! + .R^*^ < 0 

(condition (i)), are normalized as (.R^ — RW)^W = g (condition (ii”)), satisfy 

the operator inequalities > 0 (condition (iii)), and the additional conditions 

[S(®)(0)]'l'.R^®^ S*^*)(0) < [f/Wjt ^(*) ^(*) (condition (iv)) do then prove Gi{t) to be strictly 
positive, which implies Chixif) + ^min) > 0 such that —Smin is a lower bound on the 
eigenvalues of x(t)- 

As mentioned before, the vectors are not uniquely defined and become objects 
to an optimization. Good candidates are determined with a gradient-based method 
and by means of a semi-dehnite program we show that —5mm = — 10 “^ is a lower 
bound on the smallest eigenvalue of x(f) for the test case of 7 = lea, A = 2/10 0 ; and 
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/3 = 2/lOa; This HEOM induces a dynamical process with a non-Markovianity of 
A/'(Ai(t —)■ cx))) = 0.08, whose time evolution is depicted in figure 7. 

5. Conclusion 

We have demonstrated how complete positivity can be incorporated into the framework 
of HEOMs. This enables an abstraction of the equations of motion from a derivation 
based on microscopic models towards a more phenomenological perspective that opens 
a new approach to non-Markovianity. A construction of elementary models of open 
system dynamics in which physical quantities like amplitudes and frequencies of revivals 
of quantum coherence have a clear root in the underlying equations of motion - just 
like decay constants are rooted in Markovian master equations - will ultimately ease 
investigations of open system dynamics in which system size or large statistical sampling 
[31, 32, 33, 34] requires efficient phenomenological models. Although a theorem that 
describes CP dynamics as beautifully as the Lindblad theorem in the Markovian case is 
currently far out of reach for non-Markovian systems, the presented approach might be 
taken as an initial step towards such a theory. 
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Appendix A. Phenomenological constructions of HEOMs 

In section 2.4, section 3.3, and section 4.3 specific examples of HEOMs have been 
examined with respect to complete positivity. Whereas the HEOM in (40a)/(40&), 
which has been discussed in section 4.3, is a well-known HEOM for the spin-boson model 
and has first been derived in [9] based on microscopic assumptions, the other examples 
have been obtained from phenomenological considerations rather than from microscopic 
system-environment models. To this end, new techniques have been conceived that will 
be presented in Appendix A.l and Appendix A. 1.2 before the HEOMs that have been 
investigated in the main article will be derived explicitly. 

Appendix A.l. Finding HEOMs for targeted solutions 

A convenient and constructive way to obtain a well-behaved HEOM is to first specify 
the targeted system dynamics and then tailor the HEOM according to this solution. To 
this end, one first defines a dynamical map Ai(f) that propagates the system state 
Pi(f) = Ai(f) g for all times t > 0 and all initial states Pi(0) = g. The HEOM 
will then be constructed such that Ai (t) g is the first element of the extended state 


Exploring complete positivity in hierarchy equations of motion 


25 


git) = X]r=i ® K) solves the HEOM for the initial condition f?i(0) = g and 
ft(0) = O for all i >2. 

The HEOMs that we target for are of triangular form, i.e. their generators satisfy 

Cij = cjI for j = i + 1 and Cij = O for j > i + 2 (A.l) 

with \jijj being a unit of time and I and O being the identity map and the null map, 
respectively. With this form of the level k of the HEOM can affect the system 
dynamics earliest in the time-step. 

We will now successively deduce the yet unknown operators Cij, with j < i by 
imposing the condition that each operator gk{t) must vanish in the hrst k — 2 time- 
steps: Whereas the system state gi{t) does not have to vanish at all, the hrst auxiliary 
operator g 2 {t) shall vanish at the time t = 0, the second auxiliary operator g^lt) shall 
satisfy ^?3(0) = O as well as dtQ^it) |t=o = O, and this series of conditions will be continued 
for all Qkit). 

Let us demonstrate the procedure more explicitly: After having dehned a dynamical 
map Ai(t) for the system state, we dehne 

^ 2 (^ 11 , Q, t) = (j^i{t)g - £iiAi(t) g'^ /u , (A.2) 

as a function of time, the initial system state g, and the yet unknown operator £ 11 . We 
then require g 2 {Cii, g,t) to vanish at the initial time f = 0 independently of the initial 
system state g, which can be expressed by 

^^2(£ii, ^,0) = 0 Vi, j . (A.3) 

O Qij 

Solving this linear equation with respect to £11 determines the hrst level of the 
HEOM as well as the hrst auxiliary operator, which can now be obtained from 
gi(t) = £11 gi(t) -|- u g 2 (t) and reads 

02 it) = {gi{t) - £11 giit))/u . (A.4) 


In order to obtain the operators Cki to Ckk for each higher level k of the HEOM, 
we iteratively reapply this procedure. That is, we hrst dehne a function 


Qk-\-\ • • • 7 ^kk-) — I Qk(t) ^ ^ ^kj Qj (^) /' 


(A.5) 


as a function of yet unknown operators Cki to Ckk- According to the initial conditions, 
the hrst k — 2 time derivatives of ^fc+i(£fci,... ,Ckk,t) must vanish at the initial time 
t = 0, i.e. 

d ^ .X, „ . 


• • • 1 Ckki t') |t=o 0 Vi, j 


(A.6) 


is required for all n = 0 ... ,k — 2. After the operators Cki to Ckk have been determined 
such that they satisfy all these conditions, the auxiliary operator 


Qk-\-l(t^ — j Qkit^ ^ ^ j /' 


(A.7) 
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vanishes in the first k time-steps as it has been reqnired above. 

As this procednre is iteratively re-applied, the nnmber of levels of the HEOM 
increases nntil (A.7) yields Qn+i{t) = 0. This indicates that one has obtained a HEOM 
with n levels, for which the initial condition Pj(0) = O for j > 1 gives rise to the targeted 
system dynamics Qi{f) = Ai(t) pi(0). 


Appendix A. 1.1. Example: the derivation of the HEOM (21a)/(21b) The method that 
has been introdnced above is now applied in order to motivate how the HEOM in 
(21a)/(216) has been obtained. To this end, we aim for a HEOM that indnces the 
system dynamics 


Qiit) = Ai p 


1 

2 


^( 0 ) -|_ qO) jfiPj i^gO) — igiy)'j 

f{t) + ig^y'>) 


(A.8) 


with f{t) = cos(a;t) being the fnnction that characterises the reviving coherences. 
The parameters = tr(p(Ti) with i = 0,a;,|/, z and do = I determine the initial state 
g and need to satisfy = 1 and < 1 in order for g to have a 

trace eqnal to one and to be positive semi-definite. 

To find a HEOM that indnces the targeted dynamics, we parametrise the yet 
nnknown operator Cu by 


Ciig = , 


(A.9) 


where the snmmation is carried ont over = 0,x,|/, z. In accordance with (A.2), one 
defines g2{C,iiO) = (pi(^) — Qi{t))/uj snch that (A.3) reads 

0) = 0 with k = D,x,y,z (A.IO) 

and is solved with respect to Xif^- This yields = diag(—7, 0,0,7)/2 or, eqnivalently, 
-^11 2 = aT^zQI‘^ with the dephasing Lindbladian V^g = cr^gaz — g. Eqnation (A.7) 
determines the first anxiliary operator g2{t) as 


Qk{t) 


0 hkit) 

gk{t) + ig^^'>) 0 


(A.ll) 


with k = 2 and the scalar fnnction g2{t) = —e~"'^sin{ut)/2. 

We can now re-applying the same procednre to the operator g2it), in order to 

determine £21 and £22- First, these operators are parametrised in terms of yet nnknown 

scalars and xff^ as in (A.9). After defining P3(£2i,£22,= {g2{t) — JC21 gi{t) — 

£22 Q2it))/uj as in (A.5), we solve the conditions 

d d 

23(£21, £22, 0) = 0 and [^t23(£2i, £22, ^)lt=o] = ^ ’ (^-12) 

where k = 0,x,|/,z, with respect to parameters Xif^ and Xif^■ The conditions in (A. 12) 
determine only 24 ont of the 32 degrees of freedom in xfP and xff^ and the remaining 
degrees of freedom can be chosen snch that the forms of £21 and £22 become as simple 
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as possible, what eventually yields £21 Q = ojVz q /2 and £22 Q = 1 cFz QCTz- The second 
auxiliary operator Q^{t) = {Q2{t) — £21 Qi{t) — £22 Q2ii))l^ vanishes constantly such that 
the HEOM has only two level. 

We have now derived a HEOM that induces the system dynamics dehned in (A.8). 
However, we want to generalise this HEOM by permitting an arbitrary pre-factor for 
£21 such that the latter reads £21 Q = auV^ g with a free parameter a. Furthermore, 
we permit the damping constant 7 to take different values 71 and 72 on the two different 
levels of the HEOM. For a = 1/2 and 71 = 72 = 7 one obtains the original HEOM but 
the additional parameters also permit to cover different processes. The hnal HEOM 
that has now been obtained is stated in ( 21 a)/( 216 ). 

Appendix A.1.2. Example: extension of the HEOM (21a)/(21b) by means of (25) We 
can further generalise the dynamics that has been dehned in (A.8) by modifying the 
function f{t) to f{t) = -t- 2 Q![cos(a;t) — 1 ]). In the case of a = 1/2 one recovers 

the same dynamics that has been targeted for in Appendix A. 1 . 1 . However, different 
values of a permit a broader range of dynamical processes and call for a new HEOM. 

In order to obtain the latter, we follow the proposed procedure in the very same way 
that has been carried out in Appendix A. 1.1 and hnd that the hrst operator in the HEOM 
reads £11 g = 'jV^ g/ 2 . The auxiliary operator g2(t) is then determined as in (A. 7 ) with 
k = 1 and is given through (A. 11) with k = 2 and g2{t) = —e~'^*'a sin(a;t). As in (A. 5 ) we 
then dehne P3(£2i, £22, all conditions on which (compare (A.12)) are satished by the 
operators (the choice is not unique) £21 g = aw Vz g and £22 g = 5 cTz g(Jz- In contrast to 
the previous example, the operator P3(t) does in the present case not vanish constantly 
but is given through (A. 11 ) with k = 3 and g^it) = e~'^^a{l — 2 a)(l — cos(wt)). For that 
reason, an additional level in the HEOM is required, which is characterised by the yet 
unknown operators £3*, g = xff^ cn g g] with k = 1,2, 3 . By imposing the conditions 

d d d 

^24(0) = 0 , ^ [ 5 i 24 (t)li=o] = 0 > and ^ [d(h{t)\t^o\ = 0 (A. 13 ) 

with k = H,x,y,z on P4(£3i, £32, £33, t), which has been obtained from (A. 5 ) with 
k = 3 , one can eventually determine the operators £31 p = O, £32 g = w(l — 2 a) cx^ gaz, 
and £33 g = ■y^^zgo'z- Other choices for the operators are possible as (A. 13 ) does not 
determine them uniquely, but the present choice has been found most convenient. 

Eventually, the HEOM shall again be generalised by replacing the operator £32 
with £32 = (duj (Tz gcTz- The hrst two levels of hnal HEOM are then identical to ( 21 a) 
and {21b). The latter is, however, extended by an additional summand u:g^{t), which 
evolves according to ( 25 ). 

Appendix A.1.3. Example: the Jaynes-Cummings model in (18a) - (18c) In contrast to 
the previous two examples, which were not related to any particular physical system, we 
will now derive a HEOM that describes the dynamics of the resonant Jaynes-Cummings 
model [ 35 ]. To this end, we consider a two-level system inside a leaky cavity whose 
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coupling to the cavity field is determined by the coupling constant 7. The field is 
characterised by a Lorentzian spectral density that is peaked at the resonance frequency 
of the two-level system and whose spectral width is denoted by (. 

The reduced system dynamics has been obtained analytically and reads [ 35 ] 


/.(t)[g(12)]. 1 _ |/(i)Pj,(ll) • 


(A.14) 


where the parameters G M and G C determine the initial system state and 


fit) = e 


2 



+ 


— sinh 
a 



(A.15) 


with a = \/C^ — 2 'y( characterises the system dynamics. Knowing the analytical 
solution Qiit), one can apply the procedure that has been described above in order 
to find a HEOM for the resonant Jaynes-Cummings model. The derivation is carried 
out in analogy to Appendix A. 1.2 and eventually yields the HEOM ( 18 a)-( 18 c). 


Appendix A. 2 . Extending a Lindblad equation 

The HEOMs that have been tailored in Appendix A.l have been constructed under 
the premise of a previously defined solution. One can, however, also take a different 
phenomenological approach to the construction of HEOMs that is focused on the 
physical problem rather then on the mathematical form of the solution. The idea of 
this approach is to consider a Lindblad equation, which typically permits a convenient 
phenomenological interpretation, and to extend this Markovian equation to a HEOM 
that could also permit non-Markovian processes. 

As an example, we consider a two-level system that is coupled to a thermal reservoir 
with finite temperature and whose dynamics is described by the Lindblad equation 
dtQiit) = (i 3 + + BE)Qi{t) with B±Qi = 7± pi cxzp — {(j=pcr±, pi}/ 2 ) and rates 7±. 
Depending on the initial state, the dynamics that is induced by this Markovian equation 
of motion is characterised by exponential gain or loss process. 

This Markovian equation can now be extended such that it becomes a HEOM by 
adding u Q2it) to the right-hand side of the equation of motion, where l/u is a unit of 
time and Q2it) denotes the (dimensionless) auxiliary operator. The equation of motion 
for Q2if) shall be of the form dtQ2it) = C21 Pi(t) -l- C22 22(^) such that the two operators 
C21 and C22 determine to what extent the system dynamics will eventually differ from 
what would have been induced by a time-dependent Lindblad equation. How one can 
uncover an intuitive relation between the choice of £21, £22 and system dynamics is a 
challenging and in full generality yet unanswered question. A good phenomenological 
understanding can, however, be obtained for the case of £21 being proportional to £11, 
i.e. for £21 = (^/7+)(H+ + B_) with a free scalar parameter This choice of £21 causes 
the system state to affect the time-derivative dtQ2it) of the auxiliary operator in the 
same way (up to a constant factor) as the time-derivative dtQiit) of the system state 
itself. If there was no further coupling between Qi{t) and ^2(^)5 then the time-evolution 
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of ^i(t) would be completely encoded into the time-evolution of P2(t). However, due to 
the summand OJ Q2{t) in the equation of motion for pi(f), the auxiliary operator Q2{t) 
couples back to the system state and affects the latter not at the time t but only in the 
next time-step t + St. This contributes a kind of “inertia” to the system dynamics and 
gives rise to non-Markovian oscillation of quantum coherence. 

Eventually, we choose £22 such those coherences between Uaj-eigenstates and those 
between Uy-eigenstates in the auxiliary operator are damped with the rate (7+ -|-7_)/2. 
This gives rise to £22 Q = (7+ + and the HEOM that has been stated 

in ( 33 a)-( 336 ). 
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